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Extended Lagrangian Born-Oppenheimer molecular dynamics [Niklasson, Phys. Rev. Lett. 100 
123004 (2008)] has been generalized to the propagation of the electronic wavefunctions. The tech- 
nique allows highly efficient first principles molecular dynamics simulations using plane wave pseu- 
dopotential electronic structure methods that are stable and energy conserving also under incomplete 
and approximate self-consistency convergence. An implementation of the method within the plane- 
wave basis set is presented and the accuracy and efficiency is demonstrated both for semi-conductor 
and metallic materials. 
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I. INTRODUCTION 

As the available computational capacity for scientific 
computing is grovifing, first principles Born-Oppenheimer 
(BO) molecular dynamics (MD) is becoming an increas- 
ingly important tool for studying a wide range of material 
problems. First principles BOMD delivers a very accu- 
rate approach to atomistic simulations without relying 
on a fitted parameterization of the atomic interactions 
as in classical molecular dynamics. Unfortunately, appli- 
cations of BOMD that are based on self-consistent field 
(SCF) calculations such as density functional theory^"'^ 
are often limited by a very high computational cost or 
by fundamental shortcomings such as unbalanced phase 
space trajectories, numerical instabilities and a system- 
atic long-term energy drift 

Recently an extended Lagrangian BOMD (XL-BOMD) 
was introduced that avoids some of the most serious prob- 
lems of regular BOMD and enables computationally efH- 
cient and stable simulations of energy conserving (micro- 
canoncial) ensembles''"". In XL-BOMD, auxiliary elec- 
tronic degrees of freedom are included, in addition to 
the nuclear coordinates and velocities. In contrast to the 
popular extended Lagrangian Car-Parrinello molecular 
dynamics methods''^""^-, the nuclear forces are calcu- 
lated at the ground state BO potential energy surface 
and the total BO energy is a constant of motion. 

So far XL-BOMD has been limited to density ma- 
trix formulations of the extended electronic degrees of 
freedom. This excludes any practical implementation in 
widely used plane-wave pseudopotential schemes, since 
it would lead to unmanageable large density matri- 
ces. Because of the arbitrary phase of the electronic 
wavefunctions' ' it is difficult to use wavefunctions as 
the extended electronic degrees of freedom in a stable 
time-reversible or geometric integration of the equations 
of motion. Here we show how the electronic wavefunc- 
tions can be included in XL-BOMD. Our formulation al- 
lows a time-reversible integration of both the nuclear and 
the electronic degrees of freedom and it provides a highly 



efficient BOMD for plane- wave pseudopotential methods 
that is stable and energy conserving also under incom- 
plete and approximate SCF convergence. The wavefunc- 
tion XL-BOMD method was implemented in the Vienna 
Ab-initio Simulation Package (VASP)^''"^''' and its ac- 
curacy and efficiency are demonstrated both for semi- 
conductor and metallic materials. 



II. FIRST PRINCIPLE MOLECULAR 
DYNAMICS 



A. Born-Oppenheimer molecular dynamics 

First principle BOMD based on density functional the- 
ory (DFT) is given by the Lagrangian, 

£^0(R, R) = i E M,J?2 ~ [/dft[R; (1) 



where R = {Ri} are the nuclear coordinates and the dot 
denotes the time derivative. The potential C/dft [R; '^^^] 
is the ground state energy, including ion-ion repulsions, 
for the density given by the self consistent (sc) electronic 
wavefunctions, "i^^^ = {ip^j^}- Here n and k denote the 
band and reciprocal lattice vectors, respectively. The 
Euler-Lagrange equations, 
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give the equations of motion for the dynamical variables 
R(t) and R(i). 

The high cost of finding the ground state SCF solution 
'i'^^{t) is significantly reduced by using an initial guess 
that is extrapolated from previous time steps'^' 'd5'"22^ 
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In the SCF optimization in Eq. (3) above we assume a 
full optimization, which may include several iterative cy- 
cles based on, for example, simple linear mixing, Broy- 
den mixing, or the direct inversion of the iterative sub- 
space (DIIS) method" However, since the SCF opti- 
mization in practice never is complete, the extrapolation 
procedure in Eq. (3) leads to an irreversible evolution 
of the ground state electronic wave functions. The nu- 
clear forces are therefore calculated with an underlying 
electronic degrees of freedom that behave unphysically. 
The irreversibility appears most strikingly as a system- 
atic long-term energy drift By using thermostats, e.g. 
an artificial interaction with an external heat bath, these 
shortcomings of BOMD may not be noticed. However, a 
thermostat requires an underlying dynamics that is phys- 
ically correct, and the problems are therefore never re- 
moved. Only by improving the SCF convergence, which 
is increasing the computational cost, is it possible to sup- 
press the energy drift, though the problem never fully 
disappears. 



B. Wavefunction Extended-Lagrangian 
Born-Oppenheimer MD 

In our wavefunction XL-BOMD, proposed here, the 
dynamical variables of the BO Lagrangian are extended 
with a set of auxiliary wavefunctions $ = {4>nk} evolv- 
ing in harmonic oscillators centered around the self- 
consistent ground state wavefunctions '^/^'^{t), 
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Here fi is a fictitious electron mass parameter and w is a 
frequency or curvature parameter for the harmonic po- 
tentials. Applying the Euler-Lagrange equations to the 
extended Lagrangian in Eq. (4) gives 



dRi Y'dR, 



.k?dv, (5) 



motion. Equation (8) determines the dynamics of our 
auxiliary wavefunctions Since /U is set to zero, the 

only remaining undetermined parameter is the frequency 
or curvature w of the extended harmonic potentials. As 
will be shown below, lo occurs in the integration of Eq. (8) 
only as a dimensionlcss factor St^uP and therefore affects 
the dynamics in the same way as the finite integration 
time step 5t. 

Since the auxiliary wavefunctions are dynamical 
variables, they can be integrated by, for example, the 
timc-revcrsible Verlet algorithm-'. Morevover, since the 
auxiliary wavefunctions evolve in a harmonic well cen- 
tered around the ground state solution, will stay 
close 4'^'^ [t) . By maximizing the curvature of the har- 
monic extensions we can minimize their separation. Us- 
ing the auxiliary dynamical variables in the initial 
guess to the SCF optimization. 



*^'=(i) = SCF [$(t);R], 



(9) 



therefore provides an efficient SCF procedure that can 
be used within a time-reversible framework. The nu- 
clear forces will then be calculated with an underlying 
electronic degrees of freedom with the correct physical 
time-reversal symmetry. This is in contrast to conven- 
tional BOMD, where the SCF optimization is given from 
an irreversible propagation of the underlying electronic 
degrees of freedom as in Eq. (3). Hence the system will 
be propagated reversibly and should not suffer from any 
systematic drift in the total energy and phase space. 



C. Integration 

Both the nuclear and electronic degrees of freedom in 
Eqs. (7) and (8) can be integrated with the Verlet algo- 
rithm, or with other geometric integration schemes that 
preserve properties of the exact underlying fiow of the 
dynamics''*''-^. The Verlet integration of Eq. (8), includ- 
ing a weak external dissipative electronic force that re- 
moves accumulation of numerical noise'^, has the follow- 
ing form 



In the limit — > we get 

9C/dft[R; 



(6) 
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Thus, in the limit of vanishing fictitious mass parame- 
ter, ji, we recover the regular BO equations of motion 
in Eq. (7), with the total BO energy as a constant of 
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where a determines the magnitude of the dissipative force 
term with the Cm coefficients given in Ref."". The addi- 
tional electronic force introduces dissipation of numeri- 
cal noise that would accumulate in a perfectly reversible 
and lossless propagation. The dissipation breaks time- 
reversibility, but only to a high order in 5t^'^'^ . In this 
way numerical errors can be removed without causing 
any significant drift in the total energy. 
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TABLE I. Coefficients for the Verlet integration scheme with 
the external dissipative force term in Eq. (12). The coeffi- 
cients are derived in Ref.*", which contains a more complete 
set of coefficients. 
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D. Subspace alignment 

The electronic ground state wavefunctions are unique 
except with respect to their phase. This presents a prob- 
lem for the accuracy and stability of the Verlet integra- 
tion above. The wavefunctions need to be aligned to 
a common orientation to allow an accurate and stable 
integration. Aligning the wavefunctions backwards in 
time as in some previous integration schemes for regu- 
lar BOMD'^ is not possible, since it would break 
the time-reversal symmetry. We solve this problem by 
including a unitary rotation transform U in the SCF op- 
timization, which rotates ^'^'^{t) such that the deviation 
from is minimized in the Frobenius norm, i.e. 

[/ = argrmn||*"'=(t)t/' - (11) 

U can be calculated from U = (00^)~^/^0 where 
O — {'ii^'^l^) is the overlap matrix between 'ii^'^{t) and 
$(t)' Since the rotation is only applied to '^^'^{t) 
and not to previous auxiliary wavefunctions, the re- 
versibility is not affected. The redefined Verlet integra- 
tion is 

$(i + St) ^ 2$(t) - $(t - 5t) 

K 

+St^uj^{^^''{t)U - $(t)) + a ^ c,„$(t - mSt). (12) 

Note that good initial values for the axillary variables 
are important. If a poor initial guess are used the weak 
dissipation will eventually relax the auxiliary dynamics 
to a similar dynamics, but it would take time, and mean- 
while we would have bad initial guesses for the SCF opti- 
mization. In our implementation the initial values of the 
auxiliary variables, are set to the SCF optimized ground 
states, i.e. as intial conditions for $(i) we chose to set 
$(i) = ^^"^{t) for the first K + 1 time steps, where we 
perform phase alignements to the first optimized wave- 
functions ^^'^{t = to). It may be preferable to run with 
a stronger convergence criteria during the initial steps 
to get a good starting guess. In the first K + 1 initial 
steps we therefore chose to have a higher degree of SCF 
convergence than in later time steps. 



III. STABILITY AND NOISE DISSIPATION 

By aligning the phase in the SCF optimization, the 
stability of the Verlet integration in Eq. (12), under the 
condition of an approximate and incomplete SCF conver- 
gence, can be analyzed from the roots A of the charac- 
teristic equation of the homogeneous (steady state) part 
of the Verlet scheme, in the same way as for the den- 
sity matrix''''*. Assume a linearization of an approximate 
SCF optimization, Eq. (9), around the hypothetical exact 
solution where 

5f'''= = SCF[$] w -frscF(* - ^'*)- (13) 

Let 7 be the largest eigenvalue of the SCF response kernel 
FscF- Inserting Eq. (13) in the Verlet scheme, Eq. (12), 
with FscF replaced by 7, for the homogeneous steady 
state solution for which ^P* = 0, gives the characteristic 
equation 

A"+i = 2A" - A"-i + k(7 - 1)A" 

K 

+a J2 CrnA"-". (14) 

m=0 

Here the dimensionless constant k = St^oj^, and 7 S 
[— 1 , 1] is proportional to the amount of convergence in 
the SCF optimization. As long as the initial guess $(t) 
is brought closer to the ground state solution by the SCF 
procedure, I7I will be smaller than 1. If the characteristic 
roots have a magnitude |A|maa: > 1, the integration is un- 
stable (even if the accuracy is good), whereas it is stable if 
|A|maa; < 1 (eveu if the Optimization is approximate) . For 
|A|ma2: < 1 the accumulatiou of numerical noise will be 
suppressed through dissipation. By optimizing k ~ dt^uP' 
under the condition of stability under incomplete SCF 
convergence with 76 [—1,1], the curvature u"^ of the ex- 
tended harmonic wells will be maximized, which keeps 
the auxiliary wavefunctions $(i) as close as possible to 
the ground state solutions "■^^'^{t)'' . This optimization is 
performed under the additional condition of maximum 
dissipation. Our optimized values of a and k and the Cm 
coefficients can be found in Ref.^ and a few examples are 
given in Tab. I. Three different examples of dissipation 
as a function of SCF convergence as measured by |A|niax 
and I7I are shown in Fig. 1. 

IV. PLANE WAVE PSEUDO-POTENTIAL 
IMPLEMENTATION (VASP) 

Our wavefunction XL-BOMD method has been im- 
plemented in the Vienna Ab-initio Simulation Pack- 
age (VASP)^''"^'"' and the projector augmented wave 
method'*^''-. These particular methods not only require 
the integration of the wavefunctions, but also the elec- 
tron density and the Kohn-Sham eigenvalues that are 
used in the SCF optimization. In this work these addi- 
tional quantities has been be added to the Lagrangian as 



4 



a 
o 




^ 0.4 
Q 

0.2 
0.0 



BOMD, 2nd order integration 
XL-BOMD, {a=Q) 
XL-BOMD, {a>0) 



-1.0 -0.5 0.0 0.5 
SCF convergence y 



1.0 



FIG. 1. (Color online) Stability and dissipation for 2nd order 
regular BOMD' '•" and XL-BOMD, Eq. (12). The stability 
region of the regular BOMD is limited to 7 £ [—0.14,0.50] 
hence demanding a higher degree of SCF convergence, even if 
the accuracy in each step is high. In contrast, XL-BOMD is 
stable in the entire region of SCF convergence, 7 £ [—1, 1]. 



extended dynamical variables evolving in harmonic oscil- 
lators centered around their own optimized values in the 
same way as the auxiliary wavefunctions. For example, 
an auxiliary density, p(r), can be included as a dynamical 
variable through the extended Lagrangian, 

C^^°' (R, R, $, $, p, p) = C^^^iK, R, $, $) 



/ (n'^ir) - p{r)y dr. (15) 



Here p(r) follows the SCF optimized ground state density 
rf^ (r) . Also for the auxiliary density and other extended 
variables the initial values are set equal to the optimized 
ground state values for the first K + 1 steps. 
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FIG. 2. (Color online) Fluctuations in the total energy, AE, 
as a function of time for a Na bcc crystal with 16 atoms in the 
unit cell and an integration time step of 4 fs. The same SCF 
convergence criterion was used, 5E = 5 /xeV, requiering about 
2 SCF iterations per time step for both methods. The regular 
BOMD simulation shows significant systematic energy drift. 



V. APPLICATIONS 

To demonstrate the accuracy and efficiency of the 
wavefunction XL-BOMD scheme we simulate two differ- 
ent model systems with qualitatively different bonding, 
metallic sodium and semiconducting silicon. 



A. Sodium 

A unit cell of 16 bcc Na atoms was simulated for a 
total of 10,000 steps with an ionic temperature fluctuat- 
ing around 500 K with a time step of 4 fs. The regular 
BOMD integration scheme was based on a 2nd-order ex- 
trapolation of the wavefunctions from three previous time 



steps ^ 



and the XL-BOMD scheme used K = 5 ioi the 



dissipation^. Both methods used the velocity Verlet inte- 
gration for the nuclear degrees of freedom and were run 
with the same SCF energy convergence criterion, 6E = 5 
/ieV, resulting in about 2 SCF iterations per time step 
for both methods. Each SCF cycle includes one single 
construction and solution of the Hamiltonian eigenvalue 
problem. As SCF convergence accelerating algorithm we 
used the DIIS scheme-''. A plane- wave energy cutoff 
of 102 cV and a grid of 64 k-points was used and the 
exchange-correlation energy was given by the local den- 
sity approximation (LDA)'*''. 

The fluctuations in the total energy can be seen in Fig. 
2. For regular BOMD we see a small but systematic drift 
in the total energy of the order of 0.25 meV/ps. In com- 
parison, XL-BOMD shows no drift and the magnitude 
of the energy fluctuations due to the local truncation er- 
rors, occurring because of the finite time steps and the 
approximate SCF convergence, is the same. In fact, we 
have found that XL-BOMD is stable even when only 1 
SCF cycle per time step is used. This would be a general 
statement if the SCF procedure systematically improves 
the convergence in a single step ' . Unfortunately, this is 
not always the case. 



B. Silicon 

Next a Si system with 8 atoms per unit cell in a dia- 
mond structure was simulated for a total of 10,000 steps 
at an ionic temperature fiuctuating around 500 K with a 
time step of 1 fs. The SCF convergence threshold SE was 
set to 5 /ieV with a plane wave cutoff of 246 eV and a grid 
of 64 k-points was used. Otherwise the same settings as 
for sodium were applied. In Fig. 3 the fluctuations in the 
total energy AE is plotted. Also in this case, we find that 
XL-BOMD restores balance to the unphysical trajecto- 
ries of regular BOMD that shows a significant systematic 
drift in the total energy. 
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FIG. 3. (Color online) Fluctuations in the total energy, A_E 
versus time for 8 Si atoms simulated using XL-BOMD and 
regular BOMD. BOMD shows a systematic energy drift. 

VI. DISCUSSION AND SUMMARY 

In a direct comparison using the same time step and 
convergence criteria, we find that XL-BOMD and regu- 
lar BOMD have the same local truncation error, as mea- 
sured by the local amplitude of the oscillations in the 
total energy, both for the metallic and the non-metallic 
system. However, the unphysical behavior of regular 
BOMD, which has a systematic long-term energy drift, is 
removed in XL-BOMD. Only by significantly increasing 
the computational cost of regular BOMD with a higher 
degree of SCF convergence, or shorter time steps, is it 
possible to reduce the long-term energy drift. Table 
II summarizes the results of a comparison between XL- 
BOMD and regular BOMD for the Na simulation. The 
results clearly show that even though the energy drift can 
be substantially reduced in conventional BOMD, a large 
performance penalty has to be paid. XL-BOMD requires 
in general more memory, 1.5 to 2.5 times the temporary 
storage used in regular BOMD depending on the dissipa- 
tion scheme used. However, for most practical situations 
wall time is the limiting factor when running first prin- 
cipal BOMD, not memory usage. XL-BOMD therefore 
combines a more correct physical description with a lower 
computational cost. 



In many ways, XL-BOMD integrates some of the best 
features of regular BOMD and Car-Parrinello molecu- 
lar dynamics, i.e the parameter-free rigor of BOMD and 
an efficient extended Lagrangian framework as in Car- 
Parrinello molecular dynamics, where both nuclear and 
electronic degrees of freedom are included as dynamical 
variables. 

In summary, we have proposed and demonstrated a 
wavefunction XL-BOMD scheme that allows highly ef- 
ficient first principles molecular dynamics simulations 
using plane wave pseudopotential electronic structure 
methods that are stable and energy conserving also un- 
der incomplete and approximate self-consistency conver- 
gence. This extends the capability and accuracy of mod- 

TABLE II. Comparison between regular BOMD and XL- 
BOMD for a 16 atom Na bcc simulation. Time step, St, in 
fs, energy convergence threshold, SE, in /xeV, number of SCF 
cycles, and systematic Drift (per atom) in /xeV/ps. Drift < 
0.1 means no systematic drift was found. 
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ern molecular dynamics simulations. 
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